<?php

/**
 * @package JAMA
 */
require_once "../Matrix.php";

/**
 * Example of use of Matrix Class, featuring magic squares.
 */
class MagicSquareExample {

    /**
     * Generate magic square test matrix.
     * @param int n dimension of matrix
     */
    function magic($n) {

        // Odd order

        if (($n % 2) == 1) {
            $a = ($n + 1) / 2;
            $b = ($n + 1);
            for ($j = 0; $j < $n; ++$j)
                for ($i = 0; $i < $n; ++$i)
                    $M[$i][$j] = $n * (($i + $j + $a) % $n) + (($i + 2 * $j + $b) % $n) + 1;

            // Doubly Even Order
        } else if (($n % 4) == 0) {
            for ($j = 0; $j < $n; ++$j) {
                for ($i = 0; $i < $n; ++$i) {
                    if ((($i + 1) / 2) % 2 == (($j + 1) / 2) % 2)
                        $M[$i][$j] = $n * $n - $n * $i - $j;
                    else
                        $M[$i][$j] = $n * $i + $j + 1;
                }
            }

            // Singly Even Order
        } else {

            $p = $n / 2;
            $k = ($n - 2) / 4;
            $A = $this->magic($p);
            $M = array();
            for ($j = 0; $j < $p; ++$j) {
                for ($i = 0; $i < $p; ++$i) {
                    $aij = $A->get($i, $j);
                    $M[$i][$j] = $aij;
                    $M[$i][$j + $p] = $aij + 2 * $p * $p;
                    $M[$i + $p][$j] = $aij + 3 * $p * $p;
                    $M[$i + $p][$j + $p] = $aij + $p * $p;
                }
            }

            for ($i = 0; $i < $p; ++$i) {
                for ($j = 0; $j < $k; ++$j) {
                    $t = $M[$i][$j];
                    $M[$i][$j] = $M[$i + $p][$j];
                    $M[$i + $p][$j] = $t;
                }
                for ($j = $n - $k + 1; $j < $n; ++$j) {
                    $t = $M[$i][$j];
                    $M[$i][$j] = $M[$i + $p][$j];
                    $M[$i + $p][$j] = $t;
                }
            }

            $t = $M[$k][0];
            $M[$k][0] = $M[$k + $p][0];
            $M[$k + $p][0] = $t;
            $t = $M[$k][$k];
            $M[$k][$k] = $M[$k + $p][$k];
            $M[$k + $p][$k] = $t;
        }

        return new Matrix($M);
    }

    /**
     * Simple function to replicate PHP 5 behaviour
     */
    function microtime_float() {
        list($usec, $sec) = explode(" ", microtime());
        return ((float) $usec + (float) $sec);
    }

    /**
     * Tests LU, QR, SVD and symmetric Eig decompositions.
     *
     *   n       = order of magic square.
     *   trace   = diagonal sum, should be the magic sum, (n^3 + n)/2.
     *   max_eig = maximum eigenvalue of (A + A')/2, should equal trace.
     *   rank    = linear algebraic rank, should equal n if n is odd,
     *             be less than n if n is even.
     *   cond    = L_2 condition number, ratio of singular values.
     *   lu_res  = test of LU factorization, norm1(L*U-A(p,:))/(n*eps).
     *   qr_res  = test of QR factorization, norm1(Q*R-A)/(n*eps).
     */
    function main() {
        ?>
        <p>Test of Matrix Class, using magic squares.</p>
        <p>See MagicSquareExample.main() for an explanation.</p>
        <table border='1' cellspacing='0' cellpadding='4'>
            <tr>
                <th>n</th>
                <th>trace</th>
                <th>max_eig</th>
                <th>rank</th>
                <th>cond</th>
                <th>lu_res</th>
                <th>qr_res</th>
            </tr>
            <?php

            $start_time = $this->microtime_float();
            $eps = pow(2.0, -52.0);
            for ($n = 3; $n <= 6; ++$n) {
                echo "<tr>";

                echo "<td align='right'>$n</td>";

                $M = $this->magic($n);
                $t = (int) $M->trace();

                echo "<td align='right'>$t</td>";

                $O = $M->plus($M->transpose());
                $E = new EigenvalueDecomposition($O->times(0.5));
                $d = $E->getRealEigenvalues();

                echo "<td align='right'>" . $d[$n - 1] . "</td>";

                $r = $M->rank();

                echo "<td align='right'>" . $r . "</td>";

                $c = $M->cond();

                if ($c < 1 / $eps)
                    echo "<td align='right'>" . sprintf("%.3f", $c) . "</td>";
                else
                    echo "<td align='right'>Inf</td>";

                $LU = new LUDecomposition($M);
                $L = $LU->getL();
                $U = $LU->getU();
                $p = $LU->getPivot();
                // Java version: R = L.times(U).minus(M.getMatrix(p,0,n-1));
                $S = $L->times($U);
                $R = $S->minus($M->getMatrix($p, 0, $n - 1));
                $res = $R->norm1() / ($n * $eps);

                echo "<td align='right'>" . sprintf("%.3f", $res) . "</td>";

                $QR = new QRDecomposition($M);
                $Q = $QR->getQ();
                $R = $QR->getR();
                $S = $Q->times($R);
                $R = $S->minus($M);
                $res = $R->norm1() / ($n * $eps);

                echo "<td align='right'>" . sprintf("%.3f", $res) . "</td>";

                echo "</tr>";
            }
            echo "<table>";
            echo "<br />";

            $stop_time = $this->microtime_float();
            $etime = $stop_time - $start_time;

            echo "<p>Elapsed time is " . sprintf("%.4f", $etime) . " seconds.</p>";
        }

    }

    $magic = new MagicSquareExample();
    $magic->main();
    ?>
